L=8; K = 5
lambda = function(t) { 100 + 50/exp((t - 420)^2/(3600*L)) + 100/exp((L*(t - 480 - 30*L)^2)/360000)}

x = 1:1:(24*60);
# plot(x=x,y=lambda(x))

t=seq(0,24*60-1)


#2.2

# Sestrojíme K*10 hodnot z exp rozdělení pro dané lambda
event = numeric(K*10)             # Pole událostí
delta = 1 / 500                  # Delta t
tau   = 0                         # Aktuální bod na časové ose
i     = 1                         # Index

while (i<=K*10) {
  current_time = tau %/% 1        # Přiřadíme t dolní celou část tau 
  tau = tau + delta               # Tau zvýšíme o delta t
  rand = runif(1, min=0, max=1)   # Náhodná hodnota
  if (rand < lambda(current_time) * delta) {
    event[i] = tau                # S danou pravděpodobností nastane událost
    i = i + 1                     # Zvýšíme index
  }
}

eventGrid = numeric(K*10)         # Mřížka nul pro vykreslení grafu
plot (event, eventGrid)           # Vykreslení grafu časů příchodů


#2.3

per_minute = numeric(24*60)       # Vektor uschovávající počet příchodů za minutu
delta      = 1 / 1000             # Delta t
j          = 0                    # Index
tau        = 0                    # Aktuální bod na časové ose

# Bacha běží to fakt dlouho, nezvyšujte moc deltu !!!
while (j<24*60) {
  j = tau %/% 1                   # Přiřadíme j dolní celou část tau 
  tau = tau + delta               # Tau zvýšíme o delta t
  rand = runif(1, min=0, max=1)   # Náhodná hodnota
  if (rand < lambda(j) * delta) { # S danou pravděpodobností nastane událost
    per_minute[j+1] = per_minute[j+1] + 1
  }
}

write(per_minute2, file = "per_minute.txt", ncolumns = 1, append = FALSE)

# Doporučovaný histogram nepotřebuju, mám hodnoty v per_minute
plot (t,per_minute, lwd=1, type='l')   # Vykreslí experimentálně zjištěná data
lines (t,lambda(t), lwd=3, col='red')  # Přiloží graf teoretické fce
